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Abstract 



Using a discrete wavelet based space-scale decomposition (SSD), the spectrum 
of the skewness and kurtosis is developed to describe the non-Gaussian signatures 
in cosmologically interesting samples. Because the basis of the discrete wavelet is 
compactly supported, the one-point distribution of the father function coefficients 
(FFCs) taken from one realization is a good estimate of the probability distribution 
function (PDF) of density fields if the "fair sample hypothesis" holds. These FFC 
one-point distributions can also avoid the constraints of the central limit theorem on 
the detection of non-Gaussianity. Thus, the FFC one-point distributions are effective 
in detecting non-Gaussian behavior in samples such as non-Gaussian clumps embed- 
ded in a Gaussian background, regardless of the number or number density of the 
clumps. We demonstrate that the non-Gaussianity spectrum can reveal not only the 
magnitudes, but also the scales of the non-Gaussianity. 

Also calculated are the FFC one-point distributions, skewness and kurtosis spec- 
tra for real data and linearly simulated samples of QSO Lya forests. When considering 
only second and lower order of statistics, such as the number density and two-point 
correlation functions, the simulated data show the same features as the real data. 
However, the kurtosis spectra of samples given by different models are found to be 
different. On the other hand, the spectra of skewness and kurtosis for independent 
observational data sets are found to be the same. Moreover, the real data are sig- 
nificantly different from the non-Gaussianity spectrum of various possible random 
samples. Therefore, the non-Gaussian spectrum is necessary and valuable for model 
discrimination. 

Subject headings: cosmology: large-scale structure of Universe - quasars: absorption 
lines 
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1. Introduction 

This is our third paper on developing the discrete wavelet transform (DWT) 
into a powerful space-scale decomposition (SSD) method for the analysis of the large 
scale structures. The first two papers studied 1) identification of clusters (Pando & 
Fang 1996, hereafter PF1); 2) determination of power spectrum of density pertur- 
bations (Pando & Fang 1995, hereafter PF2). In this paper, a method of detecting 
non-Gaussian behavior via the DWT, and its application to both simulated and ob- 
servational data of QSO's Lya absorption lines is presented. 

There are important motivations for studying the non-Gaussianity of cosmic dis- 
tributions. In standard inflation/dark matter cosmology, primordial perturbations 
generated in the inflationary era are scale-free and Gaussian. The subsequent evo- 
lution of the density perturbations in the linear regime destroys the scaling of the 
spectrum, but the density field remains Gaussian. Deviations from the Gaussian 
state in the cosmic density distribution occur first because of the non-linear evolution 
caused by the gravitational instability of baryonic and dark matter. In this highly 
non-linear evolutionary stage the density field should be very non-Gaussian. Addi- 
tionally, even when the background density distribution is Gaussian, the distribution 
of visible objects may be non-Gaussian if their distribution is biased with respect 
to the background distribution. However, detection of non-Gaussian behavior has, in 
many respects, been inconclusive. For instance, no deviations from Gaussian behavior 
were found in the QDOT-IRAS redshift survey (Feldman, Kaiser & Peacock 1994), 
and even for simulation samples in which a strong non-linear evolution is part of the 
model, the one point distribution function was found to be consistent with a Gaussian 
distribution (Suginohara & Suto 1991). 

It is generally believed that the failure in detecting non-Gaussian behavior is 
due partially to the constraint imposed on Fourier techniques by the central limit 
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theorem when applied to the one point distribution function. It is well known from 
the central limit theorem of random fields (Adler 1981; Ivanonv & Leonenko 1989) 
that if the universe consists of a large number of dense clumps, and those clumps 
are independent, the real and imaginary components of each individual Fourier mode 
are Gaussian distributed although the probability distribution functions (PDF) of the 
field itself are highly non-Gaussian. If the clumps are not distributed independently, 
but are correlated, the central limit theorem still holds if the two-point correlation 
function of the clumps approaches zero sufficiently fast (Fan & Bardeen 1995). The 
central limit theorem holds even for processes that are a time (or space) average of 
a non-Gaussian distribution as long as the ratio of the sampling time to the scale of 
the fluctuation increases (Adler 1981). For these reasons, the one-point distribution 
function of Fourier modes is not a sensitive enough measure to detect deviations from 
Gaussian behavior. 

This difficulty can be overcome by using the count in cell (CIC) method because 
the CIC is based on localized window functions, which in essence, keep the sampling 
time to scale of fluctuation ratio mentioned above from increasing. Thus, the CIC 
statistic is not restricted by the central limit theorem. CIC has succeeded in de- 
tecting non-Gaussian signatures (Hamilton 1985; Alimi, Blanchard & Schaeffer 1990; 
Gaztanaga & Yokoyama 1993; Bouchet et al. 1993; Kofman et al. 1994; Gaztahaga 
& Frieman 1994). However, it has been found that the CIC results are dependent 
on the parameters of the window function (Juszkiewicz et al. 1993). This is because 
non-Gaussian distributions are generally scale-dependent, and therefore, a window 
function with a different scale will obtain different degrees of non-Gaussian behavior 
(Yamada & Ohkitani 1991). The results obtained via a CIC analysis will be a superpo- 
sition of information on scales larger than the size of the window. Therefore, in order 
to completely describe the non-Gaussianity of density fields and object distributions, 
it is very necessary to have an effective and uniform measure of the scale-dependence 
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or spectrum of non-Gaussianity. 

In the first part of this paper, we will study a method of detecting the non- 
Gaussian spectrum by the discrete wavelet SSD. As with the CIC, the bases of dis- 
crete wavelet SSD are localized, and therefore, avoid the restrictions of central limit 
theorem. On the other hand, the SSD modes are orthogonal and complete and no 
scale mixing occurs. It is easy to decompose the contribution of structures on various 
scales to the non-Gaussianity. Hence, one can expect that a discrete wavelet SSD will 
be effective in detecting the non-Gaussianity spectrum of large scale structures. 

In the second part of the paper, we calculate the non-Gaussianity of real and 
simulated samples of QSO Lja forests. Most 2nd order statistical studies, such as the 
two-point correlation function or power spectrum, have failed to detect structure in 
the distributions of QSO Lya forests. On the other hand, statistics not based on the 
power spectrum or 2 point correlation function have indicated that these distributions 
should have structure (see, for example, Duncan, Ostriker, & Bajtlik 1989; Liu and 
Jones 1990; Fang 1991). These examples illustrate that it is necessary to go beyond 
the 2nd order statistics, that is, to look at non-Gaussian behavior. 

2. Method 

2.1 Wavelet expansion 

We use the notation introduced in PF1 and PF2. Consider a one-dimensional 
density field p(x) over a range < x < L. It is convenient to use the density contrast 
defined by 



where p is the mean density in this field. To express S(x) in a Fourier expansion, we 
take the convention 



S(x) 



p(x) - p 
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with the coefficients computed by 

S n = I [ L 5(x)e- z2wnx / L dx. (3) 
L Jo 

Now contrast this with the wavelet expansion of the density p(x) (Meyer 1993). 
We first assume that p(x) is an L periodic function defined on space — oo < x < oo 
(this condition may actually be relaxed.). The wavelet expansion is then given by 
(Daubechies 1992, PF2) 



oo oo 



p(x)= c o , m 0o,m(a:) +pJ2 6 J,i^j,i( x ) ( 4 ) 

m=— oo j=0l=— oo 

where p is the mean density. c 0jm is the mother function coefficient (MFC) on scale 
j = and position /, while 5jj is the father function coefficient (FFC) on scale j and 
position I. The MFC and FFC are calculated by the inner products as 

/oo 
p{x)4> j>m (x)dx (5) 
-OO 

~ /"CO 

= / 5(x)ipj,i(x)dx (6) 



where 5(x) is the density contrast given by eq.(l). The mother function Jim (x) and 
father function ipjj are given, respectively, by (j>j jm (x) = (2/ ' L) l / 2 4>(2 j x/ L — m) and 
ipj,i(x) = (2 j / L) 1 ^ 2 tp(2 j x / L — I), where <f>(x) is the scaling function, and tp(x) the basic 
wavelet. The <f>(x) and ip(x) must meet certain admissibility conditions (Farge 1992). 
The Daubechies 4 wavelets meet all these conditions and are the discrete wavelet bases 
used in this paper. 

2. 2 Restriction of central limit theorem 

Among the various methods of detecting Gaussian deviations, the one-point dis- 
tribution of the density field is especially important because the probability distri- 
bution functions (PDF) of density fields can be directly determined by the one-point 
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distribution. That is, the one-point distribution can detect not only the deviation 
from a Gaussian distribution, but can also detect the non-Gaussian PDF itself. 

Let us consider the non-Gaussianity of density fields S(x) consisting of randomly 
distributed non-Gaussian clumps. In this case, eq.(3) shows that for large L the 
Fourier amplitudes, S n , are given by a superposition of a large number of non-Gaussian 
distributions. According to the central limit theorem, the distribution of 8 n will 
be Gaussian when the total number of clumps is large. Thus, in general, the non- 
Gaussianity of distributions of randomly distributed clumps can not be seen from the 
one-point distribution of the Fourier modes, 5 n , even if the PDF function of clumps 
is highly non-Gaussian. 

On the other hand, the father functions, ipj^x), are localized. If the scale of 
the clump is d, eq.(6) shows that the FFC, 8j : i, with j = \og 2 (L/d), is determined 
only by the density field in a range containing no more than one clump. That is, 
FFCs are not given by a superposition of a large number of non-Gaussian variables, 
but determined by one, or at most, a few non-Gaussian processes. The one point 
distribution of the FFC, <5^, avoids the restriction of the central limit theorem, and 
is able to detect non-Gaussianity, regardless the total number of the clumps in the 
sample being considered. 

One can study this problem from the orthonormal basis used for the expansion 
of the density field. A key condition needed for the central limit theorem to hold 
is that the modulus of the basis be less than C/y/Z, where L is the size of the 
sample and C is a constant (Ivanonv & Leonenko 1989). Obviously, all Fourier-related 
orthonormal bases satisfy this condition because the Fourier orthonormal bases in 1- 
dimension are such that (l/y/Z) | sin kx\ < C/y/Z and (l/y/Z)\ coskx\ < C/y/Z and 
C is independent of coordinates in both physical space x and scale space k. On the 
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other hand, the father functions (5) have 

liM*)l - (j) 0(1) (7) 

because the magnitude of the basic wavelet tp(x) is of the order I. The condition 
Wj,i{ x )\ < C/VT, will no longer hold for a constant C independent of scale variable 
3- 

Aside from the Daubechies 2 wavelet (i.e., the Haar wavelet), the father functions 
for the Daubechies wavelets are also localized in Fourier space. The FFC 5 j:l is only 
determined by the perturbations on scale j, regardless of perturbations on other scales. 
If the universe consists of non-Gaussian clumps on various scales j, the FFCs 5jj will 
still be effective in detecting the non-Gaussian signal. However, a measure given by 
a sum over a large number of scales will fail to do so. 

The mother functions of the discrete wavelet transform are localized in spatial 
space. But the mother functions 4>j,i{x) are not orthogonal with respect to the scale 
index j, i.e. they are not localized in Fourier space (PF1, PF2). This means that the 
MFCs are given by a sum over perturbations on all scales larger than L/2 j . Thus, 
if the number of the independent clumps in Fourier (scale) space is large, the MFCs 
will also be Gaussian by the central limit theorem. The count in cell (CIC) analysis is 
essentially the same as the MFC. The window function of the CIC corresponds to the 
mother function, and the count to the amplitude of the MFCs. Like the MFCs, the 
CIC is scale-mixed and it may not always be sensitive to the non-Gaussian behavior 
on specific scales. 

2. 3 One-point distribution of FFCs 

For the density field to be consistent with the cosmological principle, S(x) should 
be described by a homogeneous random process, i.e. its probability distribution func- 
tion should remain the same when x translated. On the other hand, the bases func- 
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tions i>j,i(x) are orthogonal to translation /. Equation (6) shows that the probability 
distribution of 5jj should be /-independent. 

The /-independence can also be seen from the relationship between FFC 5jj and 
Fourier coefficients 5 n (PF2) given as 

00 /2A~ 1/2 ~ 

5 hl = J2 7 5 n ^-n/2>)e^ (8) 

n=— 00 \ / 

where ip is the Fourier transform of ip. The wavelet basis functions ipj,i(x) are com- 
pactly supported in both x and Fourier space. Generally, ip(n) will have two symmet- 
ric peaks with centers at n = ±n p , and with width An p . The sum over integer n in 
eq.(8) need only be taken on two ranges of {n p — 0.5An p )2 j ' < n < (n p + 0.5An p )2 : ' and 
— (n p + 0.5An p )2 : ' < n < —{n p — 0.5An p )2- ? . Eq.(8) can be approximately rewritten 
as 

T s 1/2 (n p +0.5An p )2J 

1 ^ 2 

n=(n p — 0.5 An p )2J 



, j s 1/2 (n p +0.5An p )2^ 

^ ~ (-) 2 £ Re{^K)5 n e^/ 2J } 



^ v 1/2 (n p +0.5Ari p )23 

-J |^K)|2 X! \5 n \cos(e i , + e n + 2 7 rnl/2 j ) (9) 

- n=(n v — 0.5An p )2J 

where we have used ip{—n p ) = ip*(n p ) and 5_ n = 5*, because both ^(a;) and 5(x) are 
real. 6^, 6* n in eq.(9) are the phases of ip{n p ) and 5 n , respectively. As we pointed out 
in last section, 5 n is Gaussian even when the clumps are non-Gaussian. The phase 
of S n , 9 n , should be uniformly randomly distributed, and therefore, from eq.(9), the 
probability distribution of 5jj is independent of /. 

Moreover, if the spatial correlations of the random field 5 (x) decay sufficiently 
rapidly with increasing separation, ranges with different / are essentially statistically 
independent. That is, even for one realization of S(x), the values of FFCs 5jj at dif- 
ferent / can be considered as statistically independent measurements. In other words, 
each FFC can be treated as independent realizations of the stochastic variable Sjj. 
Thus, the FFCs, 5j : i, on scale j form an ensemble with 2 j realizations. The statistics 
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with respect to the one-point distribution of FFCs 5jj from one-realization should 
be equal to the results of the ensemble statistics. The goodness of this estimation is 
measured by the Large Number Theorem, that is, the relative error is about 1 / spti '. 
The one-point distribution of FFCs from one-realization will be valid in detecting sta- 
tistical features if density field S(x) is ergodic: the average over an ensemble is equal 
to the spatial average taken over one realization. 

A homogeneous Gaussian field with continuous spectrum is certainly ergodic 
(Adler 1981). More importantly, it has also been found that in some non-Gaussian 
cases, such as homogeneous and isotropic turbulence (Vanmarke, 1983), ergodicity 
also approximately holds. When one considers that the density field of the universe 
is homogeneous and isotropic, the one-point distributions of FFCs should be effective 
in measuring the non-Gaussianity of cosmic distributions S(x). 

2.4 Spectrum of non-Gaussianity 

To take advantage of the 8j/s ability to detect non-Gaussian behavior at different 
scales, j, we define the spectrum of skewness as 

1 N r 2i-l _ __ 

% = 77^3 E E [& - kifU (io) 

and the spectrum of kurtosis as 

1 N r V-l __ 
^7^EE[(^-^)1n-3, (11) 

where the variance a 2 is given by 

°i = y 3 - hi) 2 - (12) 

Eq.(12) is equal to (L/2 J )PJ ar and PJ ar used in PF2 as the spectrum of the pertur- 
bation. 

Note that eqs.(10) and (11) differ slightly from usual definition of the skewness or 
kurtosis by the sum over n. This is because at small j an individual sample covering 
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the range L will yield only a small number of 5jj. This makes the calculation of 
Kj meaningless. For instance, for j = 2, there are only two FFCs, <5 2 ,o and 5 2 ,i- In 
this case, 5jj = (62,0 + 62,1)/ 2, which using the usual definition would yield Kj = —2 
regardless of the sample or wavelet. In order to overcome this difficulty, we compile 
subsets consisting of N r samples. The number of S j: i will then be N r times larger than 
one sample making the statistics at small j viable. As with the usual definitions of 
skewness and kurtosis, Sj and Kj should vanish for a Gaussian distribution. 

3. Demonstration of non-Gaussian detection 

3. 1 Normal perturbations 

To ensure that eqs. (10) and (11) yielded the expected results for the Gaussian 
l-D density distribution p{x) was generated from Gaussian perturbations with 
the following spectrum 

P(k) = h —— (13) 

v ; 1 + 10 5 A; 4 v ; 

where k = 2irn/L, L being the range of the density field. The spectrum (13) has a 
peak at log A; ~ —1.37, or a typical scale at 1/k — 23.4 (length) units. 

Samples of distributions over L = 512 bins were produced which gave a bin size 
of 2ir units. The reconstruction of the spectrum (13) is shown in Figure la. The peak 
and the amplitude of the power spectrum are perfectly detected by the wavelet SSD 
(PF2). The results of Sj and Kj are shown in Figures lb and lc. The error bars are 



given by yJlh/N and ^96/N, respectively, where N is the total number of wavelet 
coefficients, i.e., the number of realizations times 2 j . (Press et al. 1992). As expected 
for normal perturbations, both Sj and Kj are zero. 

3.2 Distribution of clumps 

Let us consider non-Gaussian density fields consisting of clumps randomly dis- 
tribute in a white noise background. Clump distributions are often used to test 
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methods of detecting non-Gaussianity in large scale structure study (Perivolaropou- 
los 1994, Fan & Bardeen 1995). In fact, these kinds of studies have shown that one 
cannot detect the non-Gaussianity of samples by the one-point probabilities of the 
individual Fourier modes, even when the samples contain only a few independent 
clumps (Kaiser & Peacock 1991). More importantly, these distributions are necessary 
in studying whether random samples contain non-Gaussian signatures (see §3.3). 

To begin, first note that a clump or valley with density perturbation Ap c on 
length scale d at position / can be described as 

±( , = f ±Ap c if LL/2 Jc < x < (I + l)L/2 Jc 
P | otherwise 

where J c = \og 2 (L/d), and the positive sign is for a clump, the negative sign for a 
valley. If a density field p(x) consist of N randomly distributed clumps and valleys of 
scale d, so that the number density is N/2 Jc d on average, the field can be realized by 
a random variable of density perturbation 5p with a probability distribution P(Sp) 
defined as 



P[x < X] 



{ if X < -Ap c 

iV/2 Jc+1 if - A Pc < X < 

1 - iV/2 Jc if < X < A Pc 

1 if X > A Pc 



(15) 



The distribution function bp of clumps and valleys, f c (8p) can then be written ap- 
proximately as 

f c {x) = ^ = (1 - £-)5(x) + ^(x - A Pc ) + ^L5(x + Ap c ). (16) 

The 5(..) on the right hand side of eq.(16) denote 5-functions. The characteristic 

function of the random variable Sp of clumps and valleys is 

roc 2 Jc — N N 

4>c{u) = j ^ f(x)e lSu dx = -^j— + ^- cos(A cM ) (17) 

where A c = Ap c j p. It is very well known that the overall measures of skewness and 
kurtosis of the distribution (15) can be calculated from the characteristic function 
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(17). The results are 



S = 



and 



K = — 



du 3 



du 4 



(18) 



u=0 



u=0 



qJc 



3. 



(19) 



where 



a 2 = 



d 2 (j) c (u) N(A c f 

\u=Q — 



(20) 



du 2 |U=U 2 J < 
is the variance of the distribution. 

Consider density fields consisting of clumps or valleys randomly distributed in 
a background. In this case, the characteristic function is (f>(u) = (j) c (u)4>b(u), where 
4>b{u) is the characteristic function of the background distribution. For a randomly 
uniform Gaussian background with variance of, the overall variance is 



iV(A c 
2 J c 



(21) 



and the overall kurtosis is 



/2 Jc 

k= [n 



(22) 



N(s/n) 2 J ' 

where s/n = A c /cri> is the signal-to- noise ratio. Eq.(22) shows that this distribution 
becomes Gaussian when s/n is small. 

Samples of clumps and valleys randomly distributed were produced with a Gaus- 
sian noise background. Figure 2a shows a typical distribution. Here, d — 1 bin, there 
are 16 clumps and valleys on average, and the signal-to-noise ratio s/n = 5. The size 
of the distribution is L = 512 bins. The spectra, Sj and Kj, are shown in Figures 2b 
and 2c. Since the distributions are non-Gaussian, a Gaussian variance will no longer 
be applicable to estimate the statistical errors. Instead, the error bars in Figure 2 are 
calculated from the 95% confidence level from the ensemble of the samples (for each 
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parameter set at least 100 realizations are generated). The kurtosis calculated from 
the one point function of the Fourier modes has also been plotted in Figure 2. This 
figure clearly shows that non-Gaussian behavior can be detected by the FFCs, but 
not by Fourier methods. 

Figure 3 shows the kurtosis spectrum for a distribution consisting of 16 randomly 
distributed clumps and valleys with scale d = 4 bins, over a length L = 512 bins, and 
s/n — 5. Kj is significantly different from zero on scale j = 6, which corresponds 
exactly to the scale of the clumps, i.e. 512/2 6+1 = 4. (The +1 in the index is due 
to the way scaling is counted for the FFCs. See eq. (8) and figure 1 of PF1.) This 
demonstrates that the scale of the clumps can be detected by the peak in the spectrum 
of kurtosis. 

The effectiveness of detecting multiple scales of the clumps and valleys has also 
been tested. Figure 4 shows the results of generating samples consisting of 16, 32, 
and 48 clumps and valleys with a s/n = 2.0, and the scales of the clumps, d, are 
randomly distributed from 1 to 5 bins. Once again 100 realizations are generated. 
The kurtosis spectra are plotted in Figure 4. Also shown is the standard one value 
description of kurtosis, plotted for clarity of presentation at j = 9. This kurtosis is 
directly calculated from the distribution S(x) by 

^^»)-Wl-3. (23) 

Several features stand out in the Figure 4. First, the standard singled valued 
kurtosis is generally lower than that given by the wavelet kurtosis, especially, when 
the number of clumps is large. In this case the standard kurtosis totally misses 
the non-Gaussianity of the distribution. Second, the one value kurtosis contains a 
large uncertainty in detecting deviations from Gaussian behavior. This is because the 
distribution S(x) is equal to about the MFCs on finest scale. As mentioned in §2.2, 
MFC's distributions will be Gaussian if the clumps are independent on various scales. 
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On the other hand, the kurtosis spectrum detected the non-Gaussian signal at j — 6 
even when the number of clumps is as large as 48, and s/n = 2.0. Since the average 
bin width is about 3, this corresponds to about 1/3 of the 512 bins being occupied by 
clumps. The FFC's are extremely sensitive to deviations from Gaussian behavior. 

Figure 5 directly plots the one-point distributions of the FFCs for realizations with 
sixteen 4d clumps and valleys over a length 512d, and with s/n— 5. In each panel, 
the Gaussian distribution (dotted line) with the same variance and normalization as 
the FFC distribution is shown. The figure shows that the largest deviation of the 
distribution of FFCs from Gaussian distribution occurs on the scales of the clumps 
U = 6). 

3.3 Non-Gaussianity of random samples 

In large scale structure study, the usual way of generating random distributions 
covering a 1-dimensional range (xi,^) is 

Xi = xi + (x 2 - xi) ■ RAN (24) 

where Xi is the position of i-th object, and RAN is random number in (0, 1). Because 
the number of objects is an integer, the random samples given by eq.(24) easily lead 
to non-Gaussian distributions. 

In numerical calculations, the distribution S(x) is often binned into a histogram 
with a given bin size. If the bin size is less than the mean distance of neighbor 
objects, the value of the binned 5(x) will typically be or 1. The sample is then 
a d=l clump distribution with a one point distribution given by eq.(15), and not a 
Gaussian distribution. Only in the case when the mean number of objects contained 
in one bin is large does the one-point distribution approach the Gaussian case. To 
illustrate this point, Figure 6 plots the spectrum of kurtosis for a sample generated 
by eq.(24), in which the number of objects is 122 distributed in 64 bins. The figure 
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shows that the non-Gaussianity of the random sample is significant when the mean 
number of objects/bins is equal to 2. 

Similarly, the binning of real data will lead to non-Gaussianity that is obviously 
not in the original distribution. Generally, in order to maximally pick up information 
from a real data set, the bin size is taken to be the resolution of the coordinate x, 
i.e. the lowest possible size of the binning. However, this data reduction will also 
lead the lowest mean ratio of objects per bin and so consequently lead to the highest 
non-Gaussianity. 

4. Non-Gaussian detection in the Lya forests 

4-1 Simulation samples of Lya forests 

The non-Gaussian behavior of samples given by simulations of the Lya forests 
(Bi, Ge & Fang 1995, hereafter BGF) were examined. These samples have also been 
used for the demonstration of cluster identification and spectrum analysis by a wavelet 
SSD (PF1, PF2). The details of the simulations are given by BGF. The basic steps 
of the simulations are as follows: 1) generate dark matter distributions by Gaussian 
perturbations with a linear power spectrum of the standard cold dark matter model 
(SCDM), the cold plus hot dark matter model (CHDM), and the low-density flat cold 
dark matter model (LCDM); 2) generate the baryonic matter distribution by assuming 
that baryonic matter traces the dark matter distribution on scales larger than the 
Jeans length of the baryonic gas, but is smooth over structures on scales less than the 
Jeans length; 3) remove collapsed regions from the density field because Lya clouds 
are probably not virialized; 4) simulate Lya absorption spectrum as the absorption 
of neutral hydrogen in the baryonic gas, and include the effects of the observational 
instrumental point-spread-function, and along with Poisson and background noises; 
5) determine the Lya absorption line and its width from the simulated spectrum by 
the usual way of Lya line identification. 
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Within a reasonable range of the UV background radiation at high redshift, and 
the threshold of the onset of gravitational collapse of the baryonic matter, the LCDM 
model is consistent with observational features including 1) the number density of 
Lya lines and its dependencies on redshift and equivalent width; 2) the distribution 
of equivalent widths and its redshift dependence; 3) two-point correlation functions; 
and 4) the Gunn-Peterson effect. Especially important is the fact that the simulated 
data show no power in the two-point correlation function and that their 1-dimensional 
spectra is flat on scales less than 100 h _1 Mpc (BGF, PF2). 

However, no power in the two-point correlation and/or a flat spectrum does 
not mean that the sample is white noise. Instead, this may only indicate that the 
power spectrum and two-point correlation function are not suitable for describing the 
statistical features of the system being considered. Indeed, using a multi-resolution 
SSD analysis, the distributions of Lya forest samples with no power in the two- 
point correlation function have been found to be significantly different from uniformly 
random distributions on various scales (PF1). The non-Gaussianity spectrum analysis 
will support this result. 

As in PF1 and PF2, a one dimensional distribution n(x) of the Lya lines was 
formed by writing each sample into a histogram with bins on comoving scale of 2.5 
h _1 Mpc. This is about the distance at which the effect of line blending occurs. We 
then generated 100 uniformly random samples for each simulation sample via eq. (24). 
Since the lines are redshift dependent the total number of lines and the number of 
lines within a given red-shift interval (say, Az = 0.4) of the random samples were 
chosen to match the parent distribution. We calculated the FFCs and MFCs for both 
the BGF sample and the random data. The original distribution n(x) and its various 
scales of multiresoluted results can perfectly be reconstructed by the MFCs (PF1). 

Figure 7 gives the one-point distribution of FFCs for a SCDM sample (W> 
0.36A). For each scale j, the corresponding Gaussian distribution is plotted such 
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that it has the same normalization and variance as the one-point distribution. This 
figure clearly shows that all the distributions on j > 6 (or scales less than about 40 
h _1 Mpc) are significantly non-Gaussian. It is interesting to point out that the BGF 
simulation is based on a linear spectrum, and the perturbation field is Gaussian. The 
detected non-Gaussian behavior must come from the selection of high peaks in the 
Gaussian background fields (step 5 of the simulation procedure). 

Figures 8 and 9 are the spectra of skewness and kurtosis of the SCDM samples 
(W> 0.36A). In order to test whether the detected non-Gaussianity is from the 
binning, Figures 8 and 9 also show the spectra of skewness and kurtosis for two random 
samples, I and II. Random II is produced from the Gaussian one-point distributions, 
which have the same normalization and variance in each j as the SCDM sample 
(Figure 7). As expected, the kurtosis and skewness of Random I are equal to zero. 
Random I is generated using eq. (24). The redshift-dependence of the number density 
of the Lya lines is accounted for by generating the random data such that in each 
redshift range the number of lines in the random sample is the same as the SCDM 
sample. Hence, the spectra of the Random II samples is a measure of the possible 
non-Gaussianities due to the binning. 

Figure 9 shows that the amplitude of the kurtosis spectrum for the SCDM sample 
is systematically larger than the corresponding Random I sample. Recall that the error 
bars in Figure 9 do not represent the 1 a Gaussian errors, but the 95% confidence 
level from the ensemble of the samples. The difference between the spectra of the 
SCDM and Random I samples is significant. 

The skewness in Figure 8 is small, consistent with zero, but slightly positive. 
Even the Random I data has small, but positive skewness. A possible reason for the 
positive skewness is the redshift-dependence of number of Lya clouds. The FFC 
is mainly determined by the difference of (positive) densities of {j, 2k} and {j, 2k + 1} 
(PF1). Namely, for a clump in redshift space, the density change on the lower redshift 
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or lower k side contributes negative FFCs, while the higher redshift side give positive 
FFCs. If the number of Lya clumps decreases with increasing redshift, the change in 
clustering amplitudes (FFCs) on the higher redshift side (positive FFCs) should be 
less than the lower side (negative FFCs), i.e. the number of positive FFCs will be less 
than negative FFCs. In Figure 7, one can see the asymmetry (non-zero skewness) in 
the one-point distributions. 

Figure 10 and 11 give the skewness and kurtosis spectra for samples of all three 
models (W> 0.16A). For all models, the skewness is about the same however the 
kurtosis is different for different models. For the CHDM data, the Kj amplitudes are 
larger than that of the SCDM and LCDM data on all scales j. This is because there 
are far fewer high peaks in the CHDM than in SCDM and LCDM. The kurtosis is 
therefore a useful measure in distinguishing between the various models. 

4-2 Real data 

As in the first two papers (PF1 and PF2), two data sets of the Lya forests 
are examined. The first was compiled by Lu, Wolfe and Turnshek (1991, hereafter 
LWT). It contains ~ 950 lines from the spectra of 43 QSOs that exhibit neither broad 
absorption line nor metal line systems. The second is from Bechtold (1994, hereafter 
JB), which contains a total ~ 2800 lines from 78 QSO's spectra, in which 34 high 
redshift QSOs were observed at moderate resolution. In our statistics, the effect of 
proximity to z em has been considered. All lines with redshift z em > z > z em — 0.15 
were deleted from our samples. We assumed go = 1/2, so the distance of the samples 
range from a comoving distance from about D min =2,300 /i _1 Mpc to D max =3,300 
/i -1 Mpc. 

A problem in using real data to do statistics is the complex geometry of QSO's 
Lya forest. Different forest cover different spatial ranges, and no one of the forests 
distributes on the entire range of (D min , D max ). This is a difficulty in detecting the 



19 



power spectrum by some of the usual methods. At the very least, a complicated 
weighting scheme is needed. However, for the wavelet SSD, this problem can easily 
be solved. For instance, suppose a forest sample lies in a range (Di,D 2 ). The range 
can be extended to (D min , D max ) by adding zeros to the data in ranges (D min , D\) 
and (_D 2 , D max ). Since the father functions are compactly supported, the FFCs in the 
range (D±, D 2 ) will not be affected by the addition of zero in the ranges of (D m i n , D\) 
and (D 2 , D max ). Any statistic can then be computed by simply dropping all FFCs, 
iftjj, with coordinates / in the added zero ranges. The only uncertainty due to the 
boundary at (D min , D max ) are the two FFCs at boundary, ipj^ and ipj t i 2 , where l\ and 
l 2 correspond to the positions in the regions of (D m i n , D\) and (D 2 , D max ). In PF2, we 
numerically demonstrated that FFCs can correctly detect the "local" spectrum, i.e. 
spectrum in the range (Di,D 2 ), regardless of any zero padding. In other words, the 
FFC spectrum reconstruction is insensitive to the selection of boundary conditions. 
Using this technique, all samples were extended in comoving space, to cover 1024 bins 
with each bin of comoving size ~ 2.5 h _1 Mpc. Thus, all QSO samples were treated 
uniformly. 

Another problem in computing the spectrum of skewness and kurtosis of real data 
is the compilation of the subsets of the sample, needed for eqs.(10) and (11). Unlike 
the simulated samples, where as much data as needed can be generated, the available 
real data are limited. There are only N T = 43 samples (forests) for LWT and 78 for 
JB. In order to effectively use this data, M < N T files from among the complete N T 
samples are chosen to form a subset. Various combinations of the subsets M are then 
combined to form an ensemble. To investigate the effect of different combinations, 
the subsets M are formed by varying the total number of files chosen from the parent 
distribution, Nt, as well as changing the order in which the individual files are selected. 
It is found that the skewness and kurtosis calculated from these M-file ensembles are 
very stable until M contains few as 7 or 8 files, i.e. until only approximately 5% of 
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the total lines remain in the subset. The 95 % confidence intervals are then estimated 
from the ensembles. 

Figure 12 gives the FFC one-point distributions for sample of LWT > 0.36A and 
JB > 0.32A As in Figure 7, these distributions show highly non-Gaussian behavior, 
and are also asymmetric, with fewer positive FFC's and more negative ones. 

The skewness and kurtosis spectra of both the LWT and JB data are plotted 
in Figure 13 and 14, respectively. Even though the two data sets are independent 
they give the same amplitudes for the kurtosis on all scales j > 6 (or larger than 40 
h _1 Mpc). Therefore, these amplitudes are not be statistical flukes, but come from 
the real clustering of Lya clouds. Figure 14 also show the kurtosis of samples SCDM 
and LCDM. In terms of second and lower order statistics, such as number density and 
two-point correlation functions, LCDM gives best fitting of observed data. However, 
LCDM's kurtosis are found to be systematically less than real data. 

To check for non-Gaussianity due to binning, Figures 15 and 16 give the skewness 
and kurtosis spectra of the JB (W> 0.16A) sample and its Random I data. The 
Random I samples have the same number density as the Lja lines in each redshift 
ranges Az = 0.4 as the real data. The difference of kurtosis between the JB and its 
Random I sample (Figure 16) is more significant than the difference between SCDM 
and its Random I data (Figure 9). This indicates that real data is more non-Gaussian 
than the linear simulation sample. 

As in Figure 8, all the skewness of LWT, JB and JB's Random I are consistent 
with zero, but slightly positive. This is probably because the positive evolution of the 
number of Lya clusters with redshift. 

5 Discussion and conclusions 

We have demonstrated that the one-point distribution of FFCs of the discrete 
wavelet is a good tool for detecting non-Gaussian behavior in cosmic density fields. 
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The locality of the wavelet farther functions in both configuration and Fourier spaces 
allows for a way around the central limit theorem: a superposition of randomly dis- 
tributed non-Gaussian clumps will be Gaussian. 

It is generally believed that the ergodic hypothesis is reasonable if spatial cor- 
relations are decreasing sufficiently rapidly with increasing separation. In this case, 
volumes separated with distances larger than the correlation length can be considered 
as statistically independent regions. Therefore, in terms of short-range correlated 
components, such volumes can be treated as independent realizations. Many theo- 
retical models indeed predict that the perturbations in the universe are short-range 
correlated, or, at least, that the universe contains short-range correlated components. 
The wavelet FFCs effectively measure these statistically independent regions. Thus, 
in the case where the "fair sample hypothesis" (Peebles 1980) holds, the FFC one- 
point distribution taken from one-realization is a fair estimate of the PDF of density 
fields. 

The spectra of skewness and kurtosis provide a systematic and uniform measure 
of the non-Gaussianity of various samples. This method is sensitive to samples con- 
taining many clumps embedded in a Gaussian background, while Fourier methods 
fail to do so (Perivolaropoulos 1994). As opposed to the PDF given by CIC, which is 
statistically incomplete, the FFCs give a complete description of the scale-dependence 
of the skewness and kurtosis. 

The spectra of skewness and kurtosis have been detected for QSO's Lya forests in 
both observational and linearly simulated samples. In previous studies these samples 
have shown no power in the two-point correlation function and a flat power spectrum. 
The results of the non-Gaussian detection are non-trivial: all distributions are found to 
be non-Gaussian on at least scales of less than 40 h _1 Mpc, and the non-Gaussianities 
are not completely due to the effects of binning. It is clear that high order statistics, 
such as the spectrum of the kurtosis, can indeed provide information which is missed 
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by the 2nd order statistics. The amplitude and shape of the kurtosis spectrum are 
found to be the same for the two independent data sets. Thus, the features shown in 
the kurtosis spectrum should come from the formation and evolution of Lja clouds. 
We will study the dynamical implication of these features in future work. 

The kurtosis spectrum for the BGF simulation samples of Lya forests was also 
calculated. The kurtosis spectra are different for different dark matter models describ- 
ing the formation of Lya clouds. Among the BGF samples, the best fitting to the 
observed number density and its evolution of Lya clouds is given by the LCDM data. 
However, the kurtosis spectrum of the LCDM sample is significantly lower than real 
data. This result is consistent with that given by cluster identification. In PF1, it was 
found that the ratio between the numbers of larger and lower scale clusters for the real 
data is greater than that of the LCDM data. Obviously, the larger the cluster number 
ratio, the larger the deviation from a Gaussian distribution. Hence, the kurtosis and 
skewness spectrum opens a new window for looking at the statistical features of large 
scale structures. It is an important addition to the existing methods of describing the 
clustering and correlation of the cosmic density field, and for discriminating among 
models of structure formation. 

Both authors wish to thank Professor P. Carruthers, and Drs.H.G. Bi and P. Lipa 
for insightful conversations. 
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Figure captions 



Figure 1 Spectra of skewness and kurtosis for Gaussian perturbations with pertur- 
bation spectrum P(k) = k/(l + 10 5 A; 4 ), which has a peak at log A; ~ —1.37, 
where k = 2nn/L, L being the range of the sample. The samples are produced 
over L = 1024 bins, and the bin size is 2n units, a.) the reconstructed spectrum 
P{k), b.) Sj and c.) Kj. The error bars are given by ^Jlb/N for the skewness, 
and ^96/N for the kurtosis, where N is the number of wavelet coefficients. 

Figure 2 Spectra of skewness and kurtosis for clump distribution, which consist of 
sixteen d clumps and valleys randomly distributed in a Gaussian noise back- 
ground. The length of the distribution is L = 512c?, and signal-to-noise ratio 
s/n is 5. a) A typical realization of the distribution; b) Sj] c) Kj] and d) the 
kurtosis calculated from one point function of FFT. 

Figure 3 Spectra of skewness and kurtosis for clump distribution, consisting of 16 
randomly distributed clumps and valleys with scale d — 4 bins, and over a length 
L = 512 bins, s/n is 5. 

Figure 4 Kurtosis spectra for samples consisting of 16, 32, and 48 clumps. The sizes 
of the clumps, d, are randomly distributed in range of 1 to 5 bins, s/n is 2, and 
the length of the sample is L = 512 bins. 

Figure 5 Histogram of the one-point distribution of FFCs for a clump distribution 
consisting of 16 clumps and valleys with d — 4 bins in range L = 512 bins, 
and s/n = 5. The vertical coordinate is relative. At each scale j, the Gauss 
distribution (dashed line) has the same variance and normalization as the FFC 
distribution. 

Figure 6 Spectrum of kurtosis of 120 objects randomly distributed in 64 bins. 
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Figure 7 Histogram of the one-point distribution of FFCs for a BGF sample in 
SCDM model. The width of the Lya lines is > 0.36A. The vertical coordi- 
nate is the number of FFCs. At each scale j, the Gaussian distribution (dashed 
line) has the same variance and normalization as the FFC distribution. 

Figure 8 Skewness spectrum of BGF sample of SCDM model. The solid and dashed 
lines are the spectra of Random I and Random II samples, respectively. 

Figure 9 Kurtosis spectrum of BGF sample of SCDM model. The solid and dashed 
lines are the spectra of Random I and Random II samples, respectively. 

Figure 10 Skewness spectrum of BGF samples in the SCDM, LCDM and CHDM 
with W > 0.161. 

Figure 11 Kurtosis spectrum of BGF samples in the SCDM, LCDM, and CHDM 
with W > 0.16A. 

Figure 12 Histogram of one-point distribution of FFCs for samples of LWT (W> 
0.36A) and JB (W> 0.32A). The vertical coordinate is relative. At each scale 
j, the Gauss distribution (dashed line) has the same variance and normalization 
as the FFC distribution. 

Figure 13 Skewness spectrum for LWT and JB data. 

Figure 14 Kurtosis spectrum for LWT and JB data. For comparison, the corre- 
sponded Kj of the SCDM and LCDM are also plotted. The width of lines is 
> 0.16A 

Figure 15 Skewness spectra for JB data (W> 0.16A), and its Random II sample. 
Figure 16 Kurtosis spectra for JB data (W> 0.16 A), and its Random II sample. 
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